/*
 * MultiGaussian.h
 *
 *  Created on: 2011-8-17
 *      Author: yong
 */

#ifndef MULTIGAUSSIAN_H_
#define MULTIGAUSSIAN_H_
#include "../Typedef.h"
#include "MatrixUtil.h"
#include <cmath>
#if defined(_MSC_VER)
#define M_PI 3.141592653589793238462643383279502884197169399375
#endif

namespace distrim
{
template<typename T>
class MultiGuassian
{
private:
	typedef boost::numeric::ublas::matrix<T> Matrix;
public:
	MultiGuassian(const Matrix &mean, const Matrix &sigma) :
		m_dimension(mean.size2()), m_mean(mean), m_diff(1, m_mean.size2()),
				m_temp(1, m_mean.size2()), m_sclr(1, 1), m_normConst(POW(
						double(2 * M_PI), int(m_dimension / 2)) * SQRT(
						MatrixUtil::Determinant(sigma)))
	{
		MatrixUtil::LUInvertMatrix(sigma, m_inversion);
	}

	MultiGuassian(size_t dimension) :
		m_dimension(dimension), m_mean(1, dimension), m_inversion(dimension,
				dimension), m_diff(1, m_mean.size2()),
				m_temp(1, m_mean.size2()), m_sclr(1, 1)
	{
	}

	inline size_t GetDimension() const
	{
		return m_dimension;
	}

	inline void SetMeanSigma(const Matrix &mean, const Matrix &sigma)
	{
		m_mean = mean;
		MatrixUtil::LUInvertMatrix(sigma, m_inversion);
		m_normConst = POW(double(2 * M_PI), int(m_dimension / 2)) * SQRT(
				MatrixUtil::Determinant<T, T>(sigma));
	}

	T GetProbability(const Matrix &x)
	{
		m_diff = x - m_mean; // 1 * N
		prod(m_diff, m_inversion, m_temp); // 1 * N prod N * N = 1 * N
		prod(m_temp, trans(m_diff), m_sclr); // 1 * N prod N * 1 = 1 * 1.
		return EXP(-0.50 * m_sclr(0, 0)) / m_normConst;
	}

private:
	size_t m_dimension;
	T m_normConst;
	Matrix m_mean;
	Matrix m_inversion;
	// Initialize once to avoid redundant intermediate objects.
	Matrix m_diff;
	Matrix m_temp;
	Matrix m_sclr;

};
}

#endif /* MULTIGAUSSIAN_H_ */

